Mapping the global potential for marine aquaculture

Mapping global marine aquaculture potential, particularly on the West Coast of the US, focusing on optimal oyster cultivation in Exclusive Economic Zones.

Author

Affiliation

Haejin Kim

EDS223 - Geospatial

Background

Marine farming could be an important source of sustainable protein worldwide, better than land-based meat production. (Gentry et al). globally mapped potential marine farming areas, considering factors like ship traffic and oxygen levels.

Find the best Exclusive Economic Zones (EEZ) on the West Coast of the US for cultivating various oyster species. Previous research indicates that oysters thrive under specific conditions:

• sea surface temperature: 11-30°C
• depth: 0-70 meters below sea level

Datasets

Sea Surface Temperature

To characterize the average sea surface temperature in the region, use the yearly average from 2008 to 2012. The data we’re using is initially derived from NOAA’s 5km Daily Global Satellite Sea Surface Temperature Anomaly v3.1.

Bathymetry

To assess the ocean’s depth, we will utilize the following: General Bathymetric Chart of the Oceans (GEBCO).

Exclusive Economic Zones

Establish maritime borders by delineating Exclusive Economic Zones along the west coast of the United States starting from Marineregions.org.

Workflow

Below is an outline of the steps you should consider taking to achieve the assignment tasks.

library(sf)
library(dplyr)
library(spData)
library(here)
library(raster)
library(terra)
library(tmap)
library(kableExtra)

rm(list = ls())
here::i_am("index.qmd")

Clean data

To begin, load essential packages and set the path, preferably utilizing the “here” package. Proceed by reading the West Coast Exclusive Economic Zones shapefile. Next, import sea surface temperature (SST) rasters for the years 2008 to 2012, combining them into a raster stack. Additionally, read the bathymetry raster (depth.tif). Ensure that all data share a consistent coordinate reference system and reproject any datasets that deviate from the specified projection.

# load necessary package 

# set the path 
list.files(here("data"), pattern = "average *.tif", full.names = TRUE)

character(0)

# read SST rasters
sst_2008 <- rast(here( "data","average_annual_sst_2008.tif"))
sst_2009 <- rast(here("data","average_annual_sst_2009.tif"))
sst_2010 <- rast(here("data","average_annual_sst_2010.tif"))
sst_2011 <- rast(here("data","average_annual_sst_2011.tif"))
sst_2012 <- rast(here("data","average_annual_sst_2012.tif"))

# rename the column of SST rasters
names(sst_2008) <- c("temp_2008") 
names(sst_2009) <- c("temp_2009") 
names(sst_2010) <- c("temp_2010") 
names(sst_2011) <- c("temp_2011") 
names(sst_2012) <- c("temp_2012") 

# combine SST rasters 
all_sst <- c(sst_2008,
             sst_2009,
             sst_2010,
             sst_2011,
             sst_2012)

# read West Coast EEZ 
wc <- st_read(here("data", "wc_regions_clean.shp"))

Reading layer `wc_regions_clean' from data source 
  `/Users/haejinkim/Documents/UCSB/mysite/posts/2023-12-14-marine/data/wc_regions_clean.shp' 
  using driver `ESRI Shapefile'
Simple feature collection with 5 features and 5 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: -129.1635 ymin: 30.542 xmax: -117.097 ymax: 49.00031
Geodetic CRS:  WGS 84

# read in bathymetry
depth <- rast("data/depth.tif")

# check the crs
#st_crs(all_sst)  # 9122 
#st_crs(depth) # 4326 
#st_crs(wc) # 4326

# reproject using the terra way
all_sst_reproj <- project(all_sst, wc)

# check the crs
st_crs(all_sst_reproj) == st_crs(depth) # true

[1] TRUE

Process data

To proceed, we must process the Sea Surface Temperature (SST) and depth data for eventual combination. The SST and depth data possess slight variations in resolution, extents, and positions. Since we aim to maintain the integrity of the underlying depth data, we’ll employ the nearest neighbor approach to resample it and align it with the SST data.

Compute the average Sea Surface Temperature (SST) for the period 2008-2012. Convert the SST values from Kelvin to Celsius by subtracting 273.15. Adjust the extent of the depth raster to match that of the SST raster. Acknowledge the differing resolutions between the SST and depth data. Resample the depth data using the nearest neighbor approach to align with the resolution of the SST data. Confirm alignment in resolution, extent, and coordinate reference system between the depth and SST datasets. Evaluate the possibility of stacking the rasters for compatibility verification.

# compute mean SST 
mean_sst <- mean(all_sst_reproj)

# convert sst from K to C
mean_sst_c <- mean_sst - 273.15

# crop depth rast to match the extent of the SST rast
depth_cropped <- crop(depth, terra::ext(mean_sst_c)) 

# using nearest neighbor approach to resample
depth_resampled <- resample(depth_cropped,
                mean_sst_c,
                method = "near")

# stack to check if they have the same resolution
resolution_test <- c(depth_resampled,
                       all_sst_reproj) # no error message means same resolution 

Find suitable locations

To find suitable locations for marine aquaculture, reclassify SST and depth data based on oyster suitability. Set values to 1 for suitable locations and NA for unsuitable ones. Identify areas satisfying both SST and depth conditions using the lapp() function to overlay and multiply cell values.

# set suitable values 
rcl_sst <- matrix(c(-Inf, 11, NA, 
                     11, 30, 1, 
                    30, Inf, NA), 
              ncol = 3, byrow = TRUE) #anything outside of range is NA

rcl_depth <- matrix(c(-Inf, -70, NA,
                      -70, 0, 1,
                       0, Inf, NA),
                    ncol = 3, byrow = TRUE)

# reclassifying raster using a reclassification matrix
suitable_sst  <- classify(mean_sst_c,
                     rcl = rcl_sst, 
                     include.lowest = TRUE)

suitable_depth <- classify(depth_resampled,
                     rcl = rcl_depth,
                     include.lowest = TRUE)

# define function
mult_fun <- function(x, y) {
  return(x*y)}

# find locations that satisfy both conditions
overlay_suitable <- lapp(c(suitable_sst, suitable_depth), mult_fun)

plot(overlay_suitable)

Find the pixel area

This is a discovery determined manually using the Earth’s radius.

# figure out one pixel size and overall size 
res(overlay_suitable) # 0.04165905, 0.04165905

[1] 0.04165905 0.04165905

dim(overlay_suitable) # 480 x 480 

[1] 480 408   1

ext(overlay_suitable) # -131.98475233, -114.987860801091, 29.9920799888132, 49.988422964 (xmin, xmax, ymin, ymax)

SpatExtent : -131.98475233, -114.987860801091, 29.9920799888132, 49.988422964 (xmin, xmax, ymin, ymax)

# Assuming the CRS is in decimal degrees
pixel_size_degrees <- 0.04165905  # Replace this with your actual pixel size in degrees

# Convert degrees to kilometers (approximation)
latitude_degrees <- mean(c(29.9920799888132, 49.988422964))  # Mean latitude of the extent
longitude_degrees <- mean(c(-131.98475233, -114.987860801091))  # Mean longitude of the extent

# Conversion factors (approximation)
km_per_degree_latitude <- 111  # Approximation for latitude (1 degree of latitude is approximately 111 km)
km_per_degree_longitude <- 111 * cos(latitude_degrees * pi / 180)  # Approximation for longitude

# Convert pixel size from degrees to kilometers
pixel_size_km_x <- pixel_size_degrees * km_per_degree_longitude
pixel_size_km_y <- pixel_size_degrees * km_per_degree_latitude
pixel_area_km <- pixel_size_km_x *pixel_size_km_y

# Display the converted pixel sizes in kilometers
print(paste("Pixel Size in Kilometers (X):", pixel_size_km_x))

[1] "Pixel Size in Kilometers (X): 3.54281357277252"

print(paste("Pixel Size in Kilometers (Y):", pixel_size_km_y))

[1] "Pixel Size in Kilometers (Y): 4.62415455"

print(paste("Demension of Pixel Size in Kilometers (X*Y):", pixel_area_km)) # 16.382517

[1] "Demension of Pixel Size in Kilometers (X*Y): 16.3825175023378"

Determine the most suitable EEZ

We aim to assess the overall suitable area within each Exclusive Economic Zone (EEZ) to prioritize zones. To achieve this, identify suitable cells within West Coast EEZs, calculate the area of grid cells, and determine the total suitable area within each EEZ.

Find the percentage of suitability for each zone, considering rasterizing EEZ data and potentially joining the suitable area by region onto the EEZ vector data.

# making area to grid cells and mask it 
wc_rast <- rasterize(wc, overlay_suitable, 
                           field = "rgn", na.rm = TRUE)

wc_mask <- mask(wc_rast,overlay_suitable)

# find area of each suitable zones with zonal() in each of 5 regions of EEZ
## no need to use zonal 
wc_zonal <- zonal(overlay_suitable, wc_mask, na.rm = TRUE, fun = "sum")
colnames(wc_zonal)[colnames(wc_zonal) == "lyr1"] <- "count_pixel"

## use zonal
wc_zonal_summary <- wc_zonal %>% 
  group_by(rgn) %>% 
  mutate(suitable_area = count_pixel*pixel_area_km)  # pixel size converts to km^2

#join data
wc_suitable_EEZ <- full_join(wc, wc_zonal_summary , by = "rgn") %>% 
  mutate(count_pixel, 
         percentage_suitable = (suitable_area/area_km2 * 100),
         .before = geometry)

wc_suitable_EEZ %>% kable() %>% kable_minimal()

rgn	rgn_key	area_m2	rgn_id	area_km2	count_pixel	suitable_area	percentage_suitable	geometry
Oregon	OR	179994061293	1	179994.06	71	1163.1587	0.6462206	MULTIPOLYGON (((-123.4318 4...
Northern California	CA-N	164378809215	2	164378.81	11	180.2077	0.1096295	MULTIPOLYGON (((-124.2102 4...
Central California	CA-C	202738329147	3	202738.33	238	3899.0392	1.9231880	MULTIPOLYGON (((-122.9928 3...
Southern California	CA-S	206860777840	4	206860.78	197	3227.3559	1.5601585	MULTIPOLYGON (((-120.6505 3...
Washington	WA	66898309678	5	66898.31	162	2653.9678	3.9671673	MULTIPOLYGON (((-122.7675 4...

Visualize results

This visual map displays the regions where conditions are suitable for oyster cultivation. The majority of oysters thrive notably well in the southern and central regions of California. Explore these geographical patterns to gain a comprehensive understanding of oyster distribution.

#set to interactive mode
tmap_mode("view")

tmap mode set to interactive viewing

tmap_last()

Warning in tmap_last(): A map has not been created yet

NULL

#map for total suitable area for oysters by region
oysters <- tm_basemap("OpenStreetMap.Mapnik") +
  tm_shape(wc_suitable_EEZ) +
  tm_polygons(col = 'suitable_area',
              palette = 'RdYlBu',
              alpha = 0.75,
              border.col = 'black',
              title = "Total Suitable Area") +
  tm_text("rgn", size = 0.54) +
  tm_scale_bar(position = c("left", "right"))
 
oysters

This illustrates the percentage of the suitable area for oyster within each region. Conversely, Washington State exhibits a significant concentration of oyster habitats.

#set to interactive mode
tmap_mode("view")

tmap mode set to interactive viewing

tmap_last()

#map for total suitable area for oysters by region
oysters_percent <- tm_basemap("OpenStreetMap.Mapnik") +
  tm_shape(wc_suitable_EEZ) +
  tm_polygons(col = 'percentage_suitable',
              palette = 'RdYlBu',
              alpha = 0.75,
              border.col = 'black',
              title = "Total Percent of Suitable Area") +
  tm_text("rgn", size = 0.54) +
  tm_scale_bar(position = c("left", "right"))
 
oysters_percent

Peruvian Anchoveta

To enhance workflow efficiency, focus on the Peruvian Anchoveta, a species thriving within the depth range of -3 to -80 meters and temperatures spanning 13°C to 23°C. It is noteworthy that the primary habitats of Peruvian Anchoveta are located in South California, suggesting a preference for warmer waters at specific depths. Explore these precise environmental preferences to gain profound insights into the ecological distribution of Peruvian Anchoveta.

test <- function(sst_low, sst_high, depth_low, depth_high, species){

# Classify SST data: Set suitable values to 1 and unsuitable values to NA
rcl_sst<- matrix(c(-Inf, sst_low, NA, 
                     sst_low, sst_high, 1, 
                    sst_high, Inf, NA), 
              ncol = 3, byrow = TRUE) #anything outside of range is NA

rcl_depth <- matrix(c(-Inf, depth_low, NA,
                      depth_low, depth_high , 1,
                       depth_high , Inf, NA),
                    ncol = 3, byrow = TRUE)

#reclassifying raster using a reclassification matrix
suitable_sst <- classify(mean_sst_c,
                     rcl = rcl_sst, 
                     include.lowest = TRUE)

suitable_depth <- classify(depth_resampled,
                     rcl = rcl_depth,
                     include.lowest = TRUE)

#define function
mult_fun <- function(x, y) {
  return(x*y)}

#find locations that satisfy both conditions
overlay_suitable <- lapp(c(suitable_sst, suitable_depth), mult_fun)

# making area to grid cells and mask it 
wc_rast <- rasterize(wc, overlay_suitable, 
                           field = "rgn", na.rm = TRUE)

wc_mask <- mask(wc_rast,overlay_suitable)

# grid cell, extracted by wc 
extract <- terra::extract(wc_mask, wc)

# bring the pixel size 
pixel_area_km <- 16.38252

# count pixel and add up the area
extract_summary <- extract %>% 
  group_by(rgn) %>% 
  summarise(suitable_area= n()*pixel_area_km) %>% 
  filter(!is.na(rgn)) 

# join the two dataframe 
wc_suitable_EEZ <- full_join(wc, extract_summary, by = "rgn") %>% 
  mutate(percentage_suitable = (suitable_area/area_km2 * 100),
         .before = geometry) 

tmap_mode("view")
tmap_last()

#map for total suitable area by region

area_map <- tm_basemap("OpenStreetMap.Mapnik") +
  tm_shape(wc_suitable_EEZ) +
  tm_polygons(col = 'suitable_area',
              palette = 'Blues',
              alpha = 0.75,
              style = "jenks",
              border.col = 'black',
              title = paste0("Total " , species, " Suitable Area")) +
  tm_text("rgn", size = 0.54) +
  tm_scale_bar(position = c("left", "right"))

percent_area_map <- tm_basemap("OpenStreetMap.Mapnik") +
  tm_shape(wc_suitable_EEZ) +
  tm_polygons(col = 'percentage_suitable',
              palette = 'Oranges',
              alpha = 0.75,
              style = "jenks",
              border.col = 'black',
              title = paste0("Total ",species ," percent Suitable Area")) +
  tm_text("rgn", size = 0.54) +
  tm_scale_bar(position = c("left", "right"))

tmap_arrange(area_map, percent_area_map)

}

test(13, 23, -80, -3, "Peruvian Anchoveta")

tmap mode set to interactive viewing

Conclusion

This study focuses on the analysis of marine organisms beneath the ocean surface. The temperature and depth of the sea significantly impact marine life. The dataset used in this analysis provides valuable information for understanding the sea environment. While our assessments of the suitability of aquaculture are based on the present ocean conditions, it’s essential to acknowledge that the environment is currently experiencing unprecedented changes. Future efforts to evaluate how climate risks may affect aquaculture potential, considering expected shifts in regional ocean temperatures and productivity, will need to improve long-term predictions and offer more nuanced insights into how climate change will impact individual species.(Gentry et al).

Reference

1. Rebecca. R, Mapping the global potential for marine aquaculture, Nature Ecology & Evolution volume1, pages1317–1324 (2017)

Citation

BibTeX citation:

@online{kim2023,
  author = {Kim, Haejin},
  title = {Mapping the Global Potential for Marine Aquaculture},
  date = {2023-12-14},
  url = {https://khj9759.github.io/posts/2023-12-14-marine/},
  langid = {en}
}

For attribution, please cite this work as:

Kim, Haejin. 2023. “Mapping the Global Potential for Marine Aquaculture.” December 14, 2023. https://khj9759.github.io/posts/2023-12-14-marine/.